We present a new high-performance configuration interaction code optimally designed for the calculation of the lowest energy eigenstates of a few electrons in semiconductor quantum dots (also called artificial atoms) in the strong interaction regime. The implementation relies on a single-particle representation, but it is independent of the choice of the single-particle basis and, therefore, of the details of the device and configuration of external fields. Assuming no truncation of the Fock space of Slater determinants generated from the chosen single-particle basis, the code may tackle regimes where Coulomb interaction very effectively mixes many determinants. Typical strongly correlated systems lead to very large diagonalization problems; in our implementation, the secular equation is reduced to its minimal rank by exploiting the symmetry of the effective-mass interacting Hamiltonian, including square total spin. The resulting Hamiltonian is diagonalized via parallel implementation of the Lanczos algorithm. The code gives access to both wave functions and energies of first excited states. Excellent code scalability in a parallel environment is demonstrated; accuracy is tested for the case of up to eight electrons confined in a two-dimensional harmonic trap as the density is progressively diluted and correlation becomes dominant. Comparison with previous Quantum Monte Carlo simulations in the Wigner regime demonstrates power and flexibility of the method.